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ABSTRACT 


Horizontal wind components, potential temperature, and mixing 
ratio fields associated with a severe storm environment in the south 
central U.S. were objectively analyzed from synoptic upper air observa- 
tions with a non-homogeneous, anisotropic weighting function. The grid 
dimensions of the area of study were near 123 kilometers horizontally 
and 50 millibars vertically at 18 prer sure levels. Each data field was 
filtered with variational optimization analysis techniques. The vertical 
motion field was then analyzed using variational optimization analysis to 
insure that the three-component wind field satisfied mass continuity. 

The local time change of potential temperature and mixing ratio was de- 
termined in order to correlate the temperature and moisture advection 
patterns with severe storm development. A variational optimization anal- 
ysis model containing the four dimensional advection equation was used to 
produce advective forecasts of the potential temperature and mixing ratio 
fields . 

For the data base of 26/1200Z May 1973, several major findings 
emerged. The dry intrusion is characterized by warm air, the advection 
of which produces a well defined upward motion pattern, A corresponding 
downward motion pattern apparently comprising a deep vertical circulation 
in the warm air Bector of the associated low pressure system was detected. 
Another major result was the alignment of the axes of maximum dry and warm 
advection with the axis of the tornado producing squall line that subse- 
quently developed. 
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VARIATIONAL 0: 1MIZATION ANALYSIS 

OF TEMPERATURE AND MOISTURE ADVECTION IN A SEVERE STORM ENVIRONMENT 

CHAPTER I 

INTRODUCTION 

The environment of the severe thunderstorm has been subjected 
to intense study in recent years in an effort to understand the pro- 
cesses occurring within the severe thunderstorm that produce high sur- 
face winds, large hail, flood-producing heavy rains, and tornadoes. As 
the knowledge of these processes increases, so does the ability of 
meteorologists to accurately predict and perhaps eventually modify the 
occurrence and severity of these damaging phenomena. 

Empirical research, most notably by Miller (1972), provides very 
important relationships between synoptic and meso-scale features and the 
processes occurring on the individual thunderstorm scale. The basic 
reason for this is that the larger scale phenomena are continuous spa- 
tially and temporally and are subject to adequate measurement by a vari- 
ety of techniques; processes occurring on the scale of the thunderstorm 
are very difficult to measure. Since the thunderstorm's primary source 
of energy appears to be the buoyancy generated by latent heat release in 
a conditionally unstable atmosphere, studies of the moisture transport in 
a severe storm environment have been very encouraging (McGinley, 1973; 
Sasaki, 1973; Lewis, 1971; Sasaki and Lewis, 1970; and Hylton, 1972). 
These studies emphasize the high correlation between moisture convergence 
in the lower atmosphere and the occurrence and development of severe 
storms. 
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The determination of moisture, expressed as precipitable water 
(water moss per unit area) and the time change of moiBture, which can 
be correlated with precipitation (Palmen, 1967) have also been used 
with success in the analysis of thunderstorm occurrence and precipita- 
tion amounts with the storms. 

Temperature advection studies Bhould prove equally valuable. 

In addition to identifying the environmental stability changes produced 
by differential temperature advection, important clues to the dynamics 
of squall line development should result. Miller (1972) stresses the 
presence of a lower level warm air intrusion as a locator for severe 
storm development. The dynamic reasoning for this appears to be linked 
to the strong upward motions, produced by warm air advection, that are 
conducive to severe storm occurrence. Pfeffer (1962) suggests that 
thermal advection in the lower troposphere dominates over the vorticity 
advection in determining vertical motion fields for smaller scale phenom- 
ena, Morris (1972) bIbo cites the predominance of thermal advection over 
vorticity advection in producing vertical motions. 

The relationship between vertical motion, vorticity advection, 
and thickness, or thermal, advection is best described by the omega 
equation developed from the primitive equations of meteorology. The 
omega equation is a diagnostic expression for the vertical wind com- 
ponent, w, that is not reliant upon accurate measurements of the hori- 
zontal winds. Although the "tu-equation" is based on the geoB trophic 
wind approximation and the lower atmosphere has non- geos trophic com- 
ponents, an examination of ‘•he u)-equation will provide insight into the 
processes which produce vertical motion. In vector form, the ou-equation 
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simplified (from Holton, 1972) is 

“ to “ ^ [ V * 2 1) ] * C 2 r v . 7 9 ] 

where w » vertical motion, downward motion is positive 

V = horizontal vector wind 

£ “ horizontal difference operator, a vector 

Cf, ^2 “ positive constants when evaluated at constant 
pressure 

1) = absolute vortlclty, the sum of relative vorticity and 
coriolis acceleration 

9 = potential temperature 

The first term on the right hand side is the differential vor- 
ticity advection term* a term which normally predominates in the atmo- 
sphere, The second term is the thickness advection term. Upward verti- 
cal motion results from an increase with height of positive vorticity 
advection, a condition normally found from a trough axis to the down- 
stream ridge axis. This requires that the winds increase with height 
due to the ( ^ V ) tern; if the winds are invariant with height, then 
there is no contribution to vertical motion from differential vorticity 
advection. 

The thickness advection term produces upward vertical motion 
for warm air advection ( V • 7 9 is negative for warm air advection). 
This term is ordinarily small in the atmosphere because the isotherms 
parallel the wind streamlines (or height contours), but in developing 
baroclinic systems, the temperature advection term is substantial. 

The synoptic situations formed by the severe thunderstorm in the 
Great Plains all include a developing baroclinic low pressure system 
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with pronounced thermal advection patterns (Miller, 1972), Many of 
these Bystems initially develop as a warm-core dynamic low pressure 
center as a result of cyclogcnesis in the lee trough, a warm core 
trough formed from a combination of southerly warm flow, subsiding air 
across the Rocky Mountains and horizontal warm advection from the higher 
terrain to the west. The latter two mechanisms indicate that the warm 
air is very dry in the lower layers. This warm air is then advected 
eastward and northward as the low pressure center continues to develop. 
This is illustrated in Figures 1 and 2 at the 850 mb level and the 700 mb 
level for one of the investigative case studies of this paper. 

From the author's synoptic experience, the wind profile in the 
southeast sector of cyclogenesis in the lee trough is largely invariant 
with height through the lower and middle levels (occasionally to the 500 
mb level) coincident with the thermal ridge in the low levels. This sug- 
gests that the environmental vertical motion field may be largely 
determined by the thickness advection, which can be determined by eval- 
uation of temperature advection. One example of the magnitude of the 
warm air advection is shown in Figure 3, The soundingfl are for Tinker 
AFB, Oklahoma on 26/OOZ and 26/12Z, May 1973. A tornado producing squall 
line occurred in eastern Oklahoma later on in the day. 

In severe storm environments, the source for much of the heat 
energy appears to be radiational heating of the higher terrain of 
eastern New Mexico, the High Plains of TexaB and western Oklahoma. In 
severe storm synoptic situations, the cloud-free dry air to the west of 
the surface location of a maximum dew-point temperature gradient heats 
rapidly. Since the lower portions of this air mass have a near dry- 
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Figure la. 


Figure la-'lc. The 700 mb analyses for 25/1200Z, 26/1200Z, and 
27/1200Z May 1973. The synoptic patterns illus- 
trate the pronounced worm-core nature of cyclo- 
genesis and frontogenesis in the lee of the 
mountains trough, which is especially evident in 
Figure lb. During the transition from the warm 
core low to the cold core low in Figure lc, warm 
advection is pronounced to the east of the trough 
axis . This is followed by the pronounced cold 
advection to the west as the low changes to a 
cold core structure. 
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Figure lb, 



Figure lc. 
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Figure 2a. 

Figure 2a-2e. The 850 mb analyses for 25/1200Z to 27/1200Z May 

1973 at 12“hour intervals. This series illustrates 
the warm core nature of the trough in the lee of 
the mountains. As the system develops, the thermal 
ridge, which is initially coincident with the trough 
line, is advected eastward producing strong warm 
air advection. This is evident in Figures 2b to 2d . 
The source of the warm air appears to be a combina- 
tion of radiational heating from the high level 
terrain of the mountains (Figure 2a ) , southerly flow, 
subsiding air over the mountains, and remnants of the 
warm air of a subsiding ridge (not indicated in 
these figures). The strong cold air advection to 
the weBt of the developing low is associated with 
frontogeneais. 
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adiabatic lapse rate and strong winds, the vertical and horizontal 
transport of this temperature is facilitated. 

Thus, the analysis of temperature and moisture advection in 
severe storm environments is a very fruitful area of severe storm re- 
search. Studies to date have been hampered by several factors. One 
factor is the lack of measurements on the meso-scale; another is the 
accurate computation of the vertical component of motion needed for 
three dimensional analysis of the advection. It is emphasized that 
temperature and moisture advection, whether calculated or inferred from 
the synoptic situation, is a short-range forecast tool. Although ad- 
vection has a substantial role in the atmosphere (sb evidenced by the 
advection terms in the dynamic models) , it cannot describe or account 
for all the processes which occur. For example, local temperature 

(thickness) advection will produce accelerations to the wind compon- 
ents as a result of the height gradient changes. This Wind change 
would then affect the advection calculations. Due to these interact- 
ions, advection models are short range forecast models. The time 
span is limited to the order of several hours, the time scale of the 
thunderstorm and squall line. 

The advent of satellite-borne remote sensing offers, through 
profile inversion techniques based on thermal infrared measurements, 
the temperature and humidity profiles on scales even smaller than the 
meso-scale, Two residual problems are that wind measurements still 
must be derived from other sources and that measurements are contingent 
upon the availability of the satellite sensors; hence these observations 
must be time weighted before inclusion in an advection model. 
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The wind components must, at least for the near future, be 
derived from the synoptic scale upper air observation network operated 
by the National Weather Service. The measurement sites are irregularly 
spaced and the measurements include ail scales of motion and contain 
considerable noise, especially at higher levels and wind speeds (O'Brien, 
1970) . This noiBe and the non-resolvable wave lengths in the measure- 
ments must be selectively eliminated from the data. 

This research is oriented toward the development of an advec- 
tion model capable of assimilating remotely sensed temperature and 
moisture profiles (when available) with the temperature, moisture, and 
wind component data currently available from conventional observations. 
These data are objectively analyzed using an anisotropic, non-homogeneous 
weight function designed to minimize errors in the computation of the 
vertical component of motion. The data are filtered using Sasaki's 
variational optimization analysis technique to selectively suppress the 
non- resolvable wavelengths both horizontally and vertically. These 
assigned data fields are incorporated in a second variational optimiza- 
tion analysis scheme to insure that the advection quantities are contin- 
uous within the grid. The final analysis is the correlation of the 
observed fields and the linear advection of these fields with observed 
severe thunderstorm environments. This analysis model is a variational 
optimization analysis model incorporating the time dimension in addition 
to the three apace dimensions. The variational formalism includes 
filters in the four dimensions to suppress computational error and non- 
resoivable data field wavelengths introduced at times later than the 


initialization time. 
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The severe thunderstorm situation selected for study occurred 
on 26 May 1973. A prefrontal squall line developed in association with 
a developing surface low pressure center in southern Kansas. The 
squall line developed from the Kansas -Oklahoma border southward through 
eastern Oklahoma. Several tornadoes were reported with this squall 
line in Kansas and Oklahoma; one tornado in east central Oklahoma 
claimed five lives in the small town of Keefeton. From Figures 1, 2, 
and 3, the temperature and moisture advection patterns with the develop- 
ing system were extensive and pronounced; this research will correlate 
temperature and advection fields derived from variational optimization 
analysis methods with the squall line development. 

An independent data set was used to develop' and verify the nu- 
merical models, but the results are not included as a case study. This 
data set ia 11 June 1973; an areal coverage of moderate thunderstorms 
developed along the Texas -Oklahoma border. 



CHAPTER IX 


DATA OBJECTIVE ANALYSIS 

The numeric modelB of the investigation required an objective 
analysis to assign values of wind, temperature, and moisture to regu- 
larly spaced grid points from the irregularly spaced observation points. 
The objective analysis method used assigned a grid point value as the 
weighted average of the surrounding observations; the physical shape of 
the data fields determines the method of weighting the observed values. 
The weight function must be carefully selected to avoid unnecessary 
smoothing of the data fields and to avoid the introduction of frequen- 
cies resolvable by the grid network but actually non-resolvable from 
the observation point network. This may occur when the grid spacing is 
finer than the average observation point spacing. Filters may be subse- 
quently used to suppress, but not eliminate, the amplitude of non- 
resolvable wavelengths, as explaincu in Appendix B, This implies that 
the assigned data fields be as accurate as possible to avoid the intro- 
duction of error into all later numeric aspects of the investigation. 

When data to be assigned to regularly spaced grid points from 
irregularly spaced measurement points are isotropic and homogeneous, 
the weight functions may be circular and distance dependent. The radius 
of influence may be determined from the correlation function or from 
scale considerations if the analytic spectral response of the weight 
function is known. 
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However, the horizontal wind components and variables advected 
by the wind have anisotropic and non-homogeneous correlation functions. 

A correlation function, as used by Gandin (1965), is the relationship 
between the correlation coefficients of a data field and distance. 

The anisotropy stems from a greater correlation parallel to the wind 
than crosswind. The non-homogeneity arises becav.se the correlations 
are functions of the wind speed, which varies throughout the grid area. 
When an isotropic, homogeneous weight function, such as the well-known 
Cressman function, is used for the objective analysis, errors are intro- 
duced into the computations. An example of this is contained as Appendix 
A, Comparison of Cressman and Anisotropic Exponential Weighting Functions. 

An anisotropic, non-homogeneous weighting function was developed 
for this study. The weight function, W, is defined as 


w " ex P( 


-(x + v ) 


} 


k(l + P cos $) 

where V = absolute value Df the wind speed at the observation point 

V* “ defined as a scale maximum wind 
2 2 

x + y = distance squared from the observation point to the 
grid point 

k - the filter parameter 
P » the anisotropic parameter 

$ ° angle between the wind directional angle and the direc- 

tional angle from the grid point to the observation point 


This weight function was designed by combining the best features 
of Inman's anisotropic weight function (Inman, 1970) and Barnes' expo- 
nential weight function (Barnes, 1973). The exponential form is pre- 
ferred over the basic Cressman form due to the following desirable 
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features cited by Barnes j 

1. The weight function characteristics can be chosen prior to 
the analysis so that pattern scales supported by the data distribution 
will be revealed. 

2. Because the weight function approaches zero asymptotically, 
the influence of data can be extended any distance without changing the 
weight function and, therefore, the response characteristics. This 
feature, when compared to the Cressman weight function which has a 
relatively sharp cut-off at the specified maximum influence radius, is 
very desirable from a spectral response standpoint because it reduces 
the high frequencies associated with any weight function or data window. 

3. Small scale irregularities are adequately suppressed, so 
further smoothing is often not necessary. 

4. The resolution after one iteration is comparable to four 
or more iterations in the CresBman successive correction method. 

In Barnes' paper, an analytic response function of the exponen- 
tial weight function showed lc to be a low-pass filter parameter. A 
small value of k allows high frequencies in the data; a larger value of 
k suppresses these frequencies. The parameter p allows a greater weight 
to upwind /downwind observations than for crosswind observations. The 
characteristics of the weight function for selected values of k and p 
are contained in Figure 4, For the objective analysis of this investi- 
gation, k was assigned a value of 0.5 and p was assigned a value of 
4,0 after visual examination of the frequencies present in the analyzed 
wind, temperature, And moisture fields. 
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Figure 4. 


Weight assigned to station observation as s 
function of distance r, where rP = x 2 + y 2 , 
and direction vector, $, For the advection 
investigations, the grid spacing was % inch. 
Note that a family of curves is generated 
for each set of K and |3 selected, which 
governs the degree of anisotropy and scale 
smoothing. 


ORIGINAL PAGE IS 
OF POOR QUALITY 



CHAPTER III 
THE VARIATIONAL FILTER 

The assigned fields of temperature, moisture, and horizontal 
wind components require additional filtering to suppress non- resolvable 
spatial frequencies. The major requirement for filtering arises from 
the three dimensional nature of the Investigation; the initialized data 
fields have been smoothed implicitly by the weighted average of the 
observed data points in the horizontal only, and this smoothing has 
not been uniform due to the anisotropy and inhomogeneity of the weight- 
ing function. 

Since the upper air reporting network has an average spacing 
of approximately 400 km over the r ’ea of the study, the minimum resolv- 
able wavelength is 800 km under optimal conditions. In the vertical 
however, wind, temperature, and moisture data are available every 50 mb, 
the vertical minimum resolvable wavelength is 100 mb. With a computa- 
tional grid spacing of 125 km horizontally and 50 mb vertically, wave- 
lengths of less than seven horizontal grid units and two vertical grid 
units should be suppressed. 

A filter based on Sasaki's variational optimization analysis 
technique (Sasaki, 1970a, 1970b, 1971a; Wagner, 1971) 1 b used because 
the filter characteristics may be specified individually in tks? hori- 
zontal and vertical and applied simultaneously. The variational formal- 
ism is 

18 
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whert J is the functional; &, « f , and a p are the specified weak con- 
straint weights applied to the observational terms and to the first 
derivative filter terms horizontally and vertically. 

The Euler-Lagrange equations derived from equation (1) are: 

a(u-u) - + f-f ) - ® p ft = 0 (2a) 

2* ay r ap 

S(v-7) - «,< ^ ^ ) - « & - 0 (2b) 

f ax^ ay^ P dp 

These partial differential equations are elliptic and are solved in 
finite difference form by relaxation techniques. The use of variational 
optimization analysis in filtering and the derivation of the Euler- 
Lagrange equations, the solution technique, and the filter response 
characteristics are contained in Appendix B. 

Because sharp gradients on the boundary will introduce error 
into the interior grids with relaxation techniques, with a resultant 
decreased convergence rate, the six boundary faces of the three- 
dimensional grids were filtered before the filtering pass on the inte- 
rior. The boundary face filtering was accomplished with a two-dimensional 
variational filter with weak constraint weighting consistent with the 
weighting used in the interior. 

The amount of change in the fields over the 15 x 14 x 18 grid 
used, expressed as root mean square error (RMSE), for the two investi- 
gations is contained in Table 1. As shown, the filtering produced an 
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average change of less than five percent over the grid network; the 
maximum changes occurred in the vicinity of sharp gradient changes in 
the data, 

DATA FIELD 

PATE u(m sec S v (m sec S e (deg K) w(gms kg~S 

11/12Z JUN 73 .4215 .3793 .4046 .0994 

26/12Z MAY 73 .4629 .5633 .3976 .1825 

TABLE 1, Root mean square error (interior points only) for the three 
dimensional variational filtering adjustment of the two investigation 
case studies. The temperature is expressed as potential temperature, 
0, and the moisture is expressed as mixing ratio, w. 



CHAPTER XV 


VERTICAL MOTION COMPUTATION 

The variational approach will be also used to calculate the 
vertical motion and to adjust the entire wind field so that the equa- 
tion of mass continuity, necessary for advection calculations, is 
satisfied. The variational method is uniquely able to combine the 
dynamic constraint of mass continuity with the observational constraint 
to adjust the wind field such that conservation of advection quantities 
is assured. This method has been used with success by O'Brien (1970) 
and McGinley (1973), both based on Sasaki's variational optimization 
approach. The variational formalism is, in non-dimensional form, 

61=0=6 JXf {o[(u-u) 2 + (v-v) 2 ] + k(|^ + ~ + !“■) }dp dy dx (3) 

'The resulting Euler-Lagrange equations, derived by carrying out the 
indicated variations in discrete form and allowing the first variation 
to be arbitrary in the interior, are 


2^ (u-u) - - 0 

(4a) 

2a (v-v) - || = 0 

(4b) 

-£-0 

ap 

(4c) 

|h + |v + M ,o 

ox 3y op 

(4d) 
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Equations (4a) and (4b) rearranged show that the analyzed fields u and 
v are derived by adjusting the assigned fields by the gradients of the 
Lagrange multiplier, X. As a consequence of equation (4c), X is inde- 
pendent of pressure. 

The Lagrange multiplier field may be determined by the elimina- 
tion of u and v from equation (4d) by substituion of equation (4a) 
operated on by ^ and equation (4b) operated on by to produce 


2S ( a* 2 %y 2> ait 8y> 


This equation may be integrated over pressure to eliminate the terra 
containing uj, since the values of tu may be assigned at the boundaries 
of the atmosphere. Physically, the vertical motion at the top of the 
atmosphere is zero; a terrain induced vertical motion may be added at 
the base (McGinley, 1973). The Lagrange multiplier, X, is Independent 
of pressure, and the divergence of the observed wind component is 
known, so X is calculated from 


/ clA + LI 

' 2 9 

la 3x By 


-J 


( & + & ) dp 

^ ay F 


( 6 ) 


where P q = scaled pressure at base of model. This equation is an ellip- 
tic partial differential equation which can be solved for X by the re- 
laxation technique described in Appendix B. In discrete form, equations 
(5) and ( 6 ) are 

2a(v x ^ + V y v) + V~X + vjx + v p 0) = 0 (7) 

P „ 9 o 

— (7 X + v X) =2 (y u + ^v) 

2 ~ x v y p sv x v y 


( 8 ) 
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In a physical sense, the curvature of the Lagrange multiplier 

field is proportional to the deviation, or error, from mass continuity 

of the assigned wind fields. The resulting gradient of \ is then used 

to adjust the assigned fields of u and v to obtain the analyzed fields 

as shown by rearranging equations (4a) and (4b) to produce 

u = u + ~ V \ (9a) 

2a x 

V = v + ~ v \ (9b) 

2 oc y 

These analyzed u and v fields may then be used to obtain the 
analyzed vertical motion field by summation over pressure of equation 
(4d) . The vertical motion, <o, at a specified level, p, iB obtained 
through application of the trapezoidal rule to produce 

ID = u) f (v u + u v ) + (v u + y v) 1 /■o-\ 

p P+4P 2 LW x v y 'p Xv x v y 'p+Ap (9c) 

Several interesting results are obtained if mass continuity 
is treated as a weak constraint rather than a strong constraint. Al- 
though mass continuity must still be satisfied in flux calculations, 
the conversion from continuous to discrete form, the use of a vertical 
grid axis not necessarily parallel to the pressure axis (of signifi- 
cance with sloping pressure patterns), and the non-exact nature of the 
boundary conditions for omega will serve to introduce deviations from 
the calculated mass continuity. Thus, the formalism in discrete form 
is 

6J = 0 = £ 2 ff(u-u) 2 + a(v-v) 2 4- P(v x u + v v + rw) 2 (10) 

xy p y P 
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The Euler-Lagrange equations are 


o(u-u) - 0V x (v x u + V y V + “ 0 ( lla ) 

a(v-v) - Pv y (V x u + V y v + v p w) « 0 (lib) 

ev p <v + VyV + v p w) - 0 (11c) 

From equation (11c), the deviation from zero of mass-coutinuity is not 

a function of the pressure. 

If equation (11a) is operated on by Vy and equation (lib) is 
operated On by v and then subtracted from equation (Ha), the following 
result is obtained 

V y (u-u) - v x (v-v) = 0 

or, rearranged 


VV-yu = r>v-ou 
x y v x y 


( 12 ) 


or (the same result is readily obtained from equations (9a) and (9b)), 

c = e as) 


Thus, the relative vorticity, Q, of the analyzed field is identi- 
cal to the relative vorticity of the initial field, £, which implies that 
corrections made to the initial wind fields are applied only to the 
divergent component; the rotational component remains unchanged. This 
was verified in both case study investigations. With the 500 mb rela- 
tive vorticity values scaled by multiplying by 10^, the root mean square 

-5 

errors were of the order 10 . Printouts of the relative vorticity 

fields before adjustment and after adjustment at ten levels for both 
case study investigations confirmed this. This very. desirable feature 
is ordinarily assumed as a consequence of the orthogonality of the rota- 
tional and divergent components of the wind fields (O'Brien, 1970 and 
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Fankhauser, 1974). 

Solution of equations (lla-llc) is facilitated 'ey defining a 

variable 

it ° 7 u + yv + 7111 (14) 

x y p 

and rewriting the equations as 

a(u-u) - P v $ = 0 (15a) 

cf(v-v) - p VyS = 0 (15b) 

Vp$ = 0 (15c) 

The variables u and v can be eliminated from equation (14) by 

substituting equation (15a) operated on by V x and equation (15b) operated 
on by Vy to produce 

B 2 2 

$ c V + v y V + s' + Vp w ( lfi ) 

The omega term is eliminated as an unknown by summing equation 
(16) over pressure. Since $ is not a function of pressure, the analysis 
equation 

— 2 (V § + V $) - p $ = £ (V U + V v) (17) 

~ v x y o „ x v y 7 ' 

O' P ' 

is an elliptic partial differential equation, very similar to equation 
(8) (especially for very large values of R) , which is solved by relax- 
ation techniques (initial 5 is zero). The resulting field of $ is 
then used to adjust the ff and V fieldB in equations (15a) and (15b). 

The vertical motion field is derived from equation (14). 

When mass continuity is a strong constraint, the adjustment of 
the u and v fields must be accomplished through several iterations 
through the adjustment process due to the non-exact nature of the center 
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space finite difference approximations to the first derivatives used 
to evaluate the divergence and the gradients of For both case study 
Investigations, convergence of the adjusted u and v fields to the con- 
tinuity satisfied u and v fields was rapid; ten iterations were used. 



CHAPTER V 


VARIATIONAL OPTIMIZATION ADVECTION ANALYSIS 


The analysis of temperature and moisture advection in thunder- 
storm environments can be accomplished with the advection equations, 
assuming dry adiabatic processes. These equations are 


d6 

dt 
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o 
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(18b) 


Tlie temperature is expressed as potential temperature, 8, through 
Poisson's equation (Hess, 1959) 

6 <= T (1000 /p)^^ (19) 


where Rj/Cp, the ratio of the gas constant for dry air to the specific 
heat of air at constant pressure, is equal to 2/7 and p is the pressure. 
The moisture is expressed as mixing ratio, w, through the conversion of 
dew point temperature, T d , to vapor pressure, e, through Tetens empiri- 
cal formula (Saucier, 1955) 


e = 6.108 exp {17.27 T d / (237.3 + T d ) } 


(20a) 


and the conversion of vapor pressure to mixing ratio by 


w 


622 e 


(20b) 


p - e 

where •) is the pressure level, e is the vapor pressure in millibars and 
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w is the mixing ratio in grams water vapor per kilogram dry air. 

These advection equations are often applied diagnostically to 
determine advection fields, such as horizontal temperature advection 
for vertical motion analysis. Since these equations contain a time 
derivative, they are also used in prognostic form; the local time 
change is evaluated at each grid point and added to the Initial value 
to produce a forecast. A major problem with the use of the advection 
equation in prognostic form is that discontinuities appear in the fore- 


cast field, even when the wind field satisfies mass continuity. Vari- 
able propagation speeds, i.e,, the wind components, will increase the 
amount of discontinuity present. 

Computational instability is another major problem with the use 
of the advection equations as prognostic equations. When an advection 
equation is expressed in center time, center space finite difference 
approximation form in the x-t plane, the computational form iB 


%(#. 


u .At 

Xjt 


x, t+1 " ^x,t-l^ * Ax f^^x+l,t " *x~l,t^ 


( 21 ) 


where $ is the field to be advected, At is the time interval, Ax is the 
grid distance, and u is the propagation speed. As discussed in Haltiner 
(1971), this form is computationally stable for values of u At/ Ax less 
than one. When u is a variable propagation speed, the stability crite- 
rion muBt be satis field for the maximum speeds in the grid. With y and 
p grid dimension added, the stability criterion value of one for the 
coefficients decreases. 

Stability can be defined in terms of the amplification tendency 
in future time steps of a current time step error term, introduced 
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through roundoff or truncation error or from error in the data field. 
The numerical scheme is stable if the amplitude of the error remains 
the same or decreases; instability arises when the amplitude increases. 
In two dimensions, the new variable value at a time step contains 
error from only two dimensions, but with three or four dimensions, the 
error terms which may have been originally or'hogonal accumulate at 
the new time step value of the variable. Thus the stability criterion 
value decreases. For an excellent graphical presentation of the propa- 
gation of error as a function of the coefficient values, see Ketter 
and Prawel (1969) . 

Another form of instability, the computational mode, is a 

result of the use of center time finite difference approximation In 

the numerical scheme. After the first time step, normally accomplished 

with a forward time Btep (which is always unstable), the value of a 

variable, $, at the nth time step, § n , is based on advective changes 
n _2 

to $ . Thus, the values of $ at alternate time steps may diverge. 

Magata and Nishida (1971) found that this non-linear instability could 
produce the complete destruction of the advective forecast scheme, but 
could be controlled with the restart technique. They found In a simple 
advective numerical experiment that the introduction of a forward time 
step, the restart, at specified time intervals in the forecast scheme 
would control this instability. 

A final problem with the use of the advection equation as a 
forecast is that new observational data cannot be incorporated into 
the grid without producing sharp gradients which then dominate the ad- 
vective scheme. The advective forecast scheme is limited to synoptic 
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data fields; data available from remote sensing at later tinea but 
before the next synoptic observation time cannot be incorporated. 

The problems of discontinuities in the advection fields and 
the inability to update the observational data lead to the formulation 
of the temperature and moisture advection in storm environments in 
variational optimization analysis form. 

For a purely advective model, the functional is established 
as 

J « ffif C«(0-0) Z + «(w-w) 2 + o-^—) 2 + c* x <~“) 2 )dtdpdydx (22) 

The terms represent the observational constraints and the advection con- 
straints for potential temperature and mixing ratio. In this form, 
both 0 and 0 are functions of time, bo new data can be incorporated 
into the advective process at any time by adjusting the appropriate 
l) terns. 

The propagation speeds of the advection equations are the con- 
tinuity satisfied wind components to insure that the'^e are no advection 
sinkB or sources within the field. Through this requirement for mass 
continuity, the wind component fields will not vary with time at the 
grid points (this can be established in a functional with observational 
constraints on the wind components and with the local time change 
included as either a weak or strong constraint). Ihe wind components 
remain variable with respect to the spatial derivatives and thus take 
part in the variations with interesting consequences * 

For simplicity, consider a functional similar to equation (22) 
but for 0 in the x-t dimensions only. 
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J - JT (S(0-8) Z + dt dx (23) 

6J ° o ■ 6 Jfee-e) 2 + ^<|f + u dx ( 24 > 

0 - jj(2'3(e-'e)60 + 20 1 (|- + u I”) 6(|| + U ||)}dt dx (25) 

- J J{ 2 cy<e-'e > 66 + 2 of 1 C 60 (- -j| - |j(u |f)) ] 

+ 2a l 6u ff dt 3 dt dx (26) 


The Euler-Lagrange equation from the variation of u is 




(27) 


Since ^ is non-zero, then d0/dt is exactly zero, which is analogous to 

OX 

the use of the advcetion equation as a strong constraint; the local 
time change of 0 is determined solely by advective changes. 

The derivation of the Euler-La grange equations from equation 
(22) written for 0 is 

6J » 0 = 6 JJJWV +a l <$) 2 } dt dp dy dx (28) 

* ifjjc 2 «(0-0) 6e + 2a i if ( ^ft + 6u n + u 


+ 6v ^ + v + 6ou ^ + u) ) dt dp dy dx 
dy dy dp ap 


, d0 v 

at dt ax dt 


e JJlf[6e£(e-e) - «i ^ 

-£< V SS> + 


(29) 


+ 6vcv i If df 6ujq, i ft 5 dt dp dy dx 


(30) 
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The ocher problem is the difficulty of numerical solution of 
this parabolic equation with four dimensions. With the time dimension 
plus any other spatial dimension, a marching method, common to the 
numerical solution of parabolic equations, will generate three unknown 
coefficients at each time step due to the finite difference approxi- 
mation form of the mixed derivative term. At each time step, a square 
tridiagonal matrix of these unknown coefficients is generated. When 
the computational stability criteria are satisfied, the tridiagonal 
matrix is well- conditioned (dominant diagonal elements) and the coeffi- 
cients are easily determined by Gauss elimination (Ketter and Prawel, 
1969). When other spatial dimensions are added, the numerical solution 
becomes very complex. 

Both of these problem areas can be adequately resolved by the 
addition of temporal and spatial first derivative filters to the func- 
tional (Sasaki, 1971b) to produce (in 8 only for simplicity) 

j = UJJ { ate-e ) 2 + «,< jff > 2 + Si C |f > 2 

+ ff > 2 + ( If > 2 ) + %< H > 2 }dt; dp dy dx (34) 

The inclusion of the filter terms effectively suppresses the 
discontinuities, which are high frequency. The solution field of 8 is 
basically an advection field with the discontinuities suppressed. As 
a very essential additional benefit, the inclusion of the filter terms 
changes the analysis equations from parabolic partial differential 
equations to the more tractable elliptic form. The Euler-Lagrange 
equation from the variation of 8 is now 

iiiECEDiNq page! blank Not pilmedi 
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c^C^-f) + 2-2 + p 2 (--jr + ^-|) + P 3 (^-f) a a(0-9) (35) 

1 dt 3t Z Z 3x' 3y J 3p 

Equation (35) expanded and rearranged forma the analysis equa- 


3 + <“i“ 2 + ?2> 3 + (V + 
at ax ay 


+ fajn? tp 3 )^|n«i < u !dh + v Ip 

^ + uv|?|- + u«|^ + v.fc|-) 


s!0- +v + 
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. „ r „^Mi av 30 , Suj M \ a.,, / du ae 
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By 3y ay dp' ' an ^x rto an an an ' J 


dp ax a? dp ap ap 


= a (0 - B) (36) 

This equation is an elliptic type partial differential equation, which 

Is solved through relaxation techniques. The type of equation may be 

verified in any two dimensions by examining the characteristic quadratic 

2 

form (the familiar B -4AC < 0) . In x and t. 


4u 2 c^ -4 (a l + P t ) (0 2 + C^u 2 ) < 0 

is satisfied for any positive weights o^, 0^, and 0 2< Note that without 

the addition of the weights (the filter terms), equation (36) reduces 

to a parabolic type equation. 

. 2 2 .22. 

4u oc^ - 4u = 0 

Equation (35), a linear partial differential equation with vari- 
able coefficients, is the primary analysis equation for the advection 
of temperature; a comparable equation for moisture is obtained by the 



same derivation technique. 


The solution of equation (35) requires that the "observed" 
fields of temperature and mixing ratio be established for all time 
steps. As a result of the variational optimization analysis, which 
determines the best fit of a field subject to the constraints imposed, 
the observed data fields must be as accurate as possible. Accuracy in 
this sense implies a consistency of the observed fields with the diag- 
nostic and prognostic constraints, although a simultaneous consistency 
is not possible without using a variational model. T.jis indicates that 
the observed data fields in time should be determined by advection from 
the initial time. The result should be a faster convergence (savings 
on computer time) to more realistic final data fields. 

An interesting corollary to this point is that perhaps the most 
Important application of the variational optimization analysis method 
is to adjust observed data fields to a form consistent with the diag- 
nostic and prognostic equations of the model. The data fields are 
adjusted using explicitly all available information, which may be sta- 
tistical, dynamic, kinematic, thermodynamic, or energetic (see, for 
example, Lewis and Grayson, 1972). 

The observed fields at the time steps were determined by the 
linear advection equations in the form 
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(37) 


In scaled, non-dimensional form using center time and center space 
finite difference approximations, equation (37) is 

- V 9 - (u'V 0 + v’V 9) + oj'V 6 

t Ax x y Ap p 


( 38 ) 
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where 


V* = 10 ra sec 


-1 


W* - 10 ^ mb sec ^ 


u = u'V* , 
v = v'V* 

<u = cu'W* , 

5 

A* « horizontal grid distance, approximately 1,25 x 10 m 
Ap = vertical grid distance, 50 mb 
At *= time grid scale, 600 sec 

n n+1 n-1' 


In view of the variable propagation speeds, the scheme was re- 
started every two hours (7200 sec) as a precaution to prevent the high 
frequency computational instability. 



CHAPTER VI 


RESULTS 


With cyclogenesis and frontogenesis in the lee-of-the-mountaits 
trough, temperature and moisture advection patterns from the surface 
to the 700 millibar level are pronounced, as shown in Figures 1, 2, 
and 3, From the omega equation, upward motion should be associated 
with the low level temperature advection and from Miller (1972) , storm 
development should be associated with the moist and dry advection pat- 
terns. He states, 

In situations preceding significant tornado 
development, a distinct dry tongue is present 
in low or middle levels, and, provided other 
criteria are satisfied, the primary develop- 
ment will occur where the dry tongue intrudes 
into or over the lower moist tongue... 

Dry air intrusions not only help in delineating 
future tornado and severe weather areas, but 
apparently provide a major contribution to the 
trigger mechanism in the majority of tornado 
situations . 


The 26/1200Z May 1973 case study shows that the warm air of the 
thermal ridge is also very dry; most dew point depressions in the dry 
air were reported as the maximum reportable 30 degrees Celsius, Several 
figures emphasize these features and that they are of sufficient hori- 
zontal and vertical extent to be resolvable in the advection analysis. 
Figure 5a shows the pressure and frontal patterns at 26/1200Z May 1973 
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Figure 5a. NMC Surface Analysis, 26/1200Z MAY 1973 



Figure 5b. NMC Surface Analysis, 26/2100Z MAY 1973, and 
Radar Summary, 26/2040Z MAY 1973. The shaded 
area indicates thunderstorm line activity; 
cells with tops over 40,000 ft. are indicated 
by large dots. 
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and Figure 5b shows the pressure and frontal pattemB at 26/2100Z with 
the radar detected thunderstorm coverage at 26/2040Z. At 26/1200Z, 
there were no radar echoes in Oklahoma except along the extreme north- 
ern border with Kansas. The prefrontal squall line depicted in Figure 
5b originated in north central Oklahoma and developed southward. The 
radar echo coverage included for reference in succeeding figures is 
the 26/2140Z thunderstorm activity. 

As shown in Figure 3, the warm, dry air lies over cooler, 
moist air with a very sharp inversion surface at the boundary. Figure 
6 shows that this inversion is horizontal over several hundred kilo- 
meters extent; the location where this sharp moisture discontinuity 
intersects the higher terrain of West Texas is known as the dry line, 
dry front, or Marfa front (from Marfa in southwest Texas). The dry 
line location at 26/1200Z was between Midland TX and Abilene TX; the 
Midland dew point temperature was -5 deg C (23 deg F) and the Abilene 
dew point was 21 deg C (70 deg F) . The temperature and dew point pro- 
file for these two stations is shown in Figure 7; the moist layer at 
Abilene is shallow. The high dew point temperature at Abilene is part 
of the surface layer moisture ridge, which is shown in Figure 8. 

The advection analysis grid was a 11 by 10 by 18; the grid 
spacing was near 125 km horizontally and 50 mb vertically. Twelve 
time steps were used in the four dimensional advection analysis; the 
time interval was 20 minutes. The actual wind speeds were used as the 
advective propagation speeds. The horizontal grid network is shown in 
Figure 9. 

The continuity satisfied vertical motion field at 26/12Z 
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Figure 6. Temperature profiles for Abilene TX (ABI) , 
Greater Southwest Airport TX (GSW), and 
Tinker AFB OK (TIK) 26/1200Z May 1973. 
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Figure 7 » Temperature and dew point temperature 

profiles for Abilene TX (ABI) and Midland 
TX (MAF), 26/1200Z May 1973. The surface 
dew point at Midland is -5 deg C (23 deg F) 
and at Abilene is 21 deg C (70 deg F) . 
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Figure 8. The dew point temperature at 26/1200Z 
May 1973. The sharp gradient in 
western Texas represents the dry line. 
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Figure 9« The 11 x 10 horizontal grid established for 
the advection analysis. The vertical grid 
has 18 levels at 50 mb intervals; the hori- 
zontal grid spacing is near 125 km. Section 
A-A' is the cross section used in the 26/1200Z 
MAY 1973 analysis. 
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contained several features of interest. The omega fields at the 900, 
800, 700, 600, and 500 millibar levels are in Figure 10. The axis of 
maximum upward vertical motions (negative w) coincides with the area 
of maximum temperature advection apparent in Figure lb and 2c. This 
axis also bears a strong similarity to the squall line location. A 
secondary axis of upward vertical motion is parallel to the frontal 
thunderstorm development. This is rooBt evident at the 700 millibar 
level. At the 500 millibar level, an area of downward vertical motion 
appears through central Oklahoma, Downward motion behind a squall line 
has been proposed as a compensatory motion which transfers the higher 
wiml speeds to lower levels (Miller, 1972) ; the appearance of this 
downward motion may be purely coincidental, although it does occur in 
these data before squall line initiation. 

The vertical motion in the cross-section along line A-A’ (see 
Figure 9) perpendicular to the squall line axis shows that the downward 
vertical motion at the 500 millibar level iB part of an organized circu 
lation, illustrated in Figure 11. The nature of this potentially very 
significant circulation pattern will have to be verified in other case 
studies; perhaps the downward motion is a compensating motion for the 
upward motion associated with the warm air advection at lower levels. 
This would produce warming and drying with the subsiding air, which 
could then enhance the warm, dry intrusion. 

The existing temperature and moisture advection fields were cal 
culated at 26/1200Z from the advection equations 
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Figure 10. Vertical motion, w x 10 mb/sec, at 26/1200Z 
May 1973. The radar sundry at 26/2140Z May 
GE IS ^973 is Included for reference. 






DI'TAHCt ItCTIOW A-A' 


3 

Figure 11. Vertical motion, u x 10 mb/sec, along 
the cross section A-A'. The downward 
motion areas (positive ui) are Indicated 
by the symbol +. 
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These udvective fields, Figure 12 and 13, show that the maximum tempera- 
ture advection is positively correlated with the maximum drying. The axes 
of the maximum drying and maximum warming are both parallel to the squall 
line, as anticipated by application of Miller's (1972) empirical rules. 

Figure 13 also illustrates the low level moisture increase at the 900 
millibar level beneath the drying associated with the warm, dry intrusion. 
This pattern is more pronounced in cross-section in Figure 14, which includes 
the 950 millibar level advective fields. The moisture increase at low level 
on the left hand side of the section is associated with the moisture advec- 
tion from the moisture ridge depicted in Figure 8. The low level tempera- 
ture advection fields are weak bfeneath the warm intrusion, as shown in 
cross section in Figure 15. 

The advective fields computed from the forecast potential tempera- 
ture and mixing ratio fields show the continued development of the patterns 
Miller cites as necessary for severe storm development, as shown in Figures 
16, 17, 18, and 19. The advective patterns reflect moisture advection be- 
low the warm, dry intrusion, which is being advected eastward. This is 
coincident with continued weak moist advection above the dry intrusion, 
not shown in these figures but evident from Figure 3. This pattern is 
cited by Miller as a possible trigger mechanism for the squall line initia- 
tion, as shown in Figure 20. The temperature advection pattern shows the 
strengthening of the mid-level cooling behind the intrusion and continued 
warm advection with the intrusion. 

As suggested by Miller and as indicated by the results of the advec- 
tion model (Figures 12 and 16), the severe thunderstorms initially develop 
through the warm, dry inversion layer. Further evidence of this is contained 
in the ATS-III imagery of 26 May 1973. The ATS III imagery and possible 
interactions of the mesoscale warm, dry intrusion and the thunderstorm scale 
processes are discussed in Appendix C. 
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Figure 12, The local change of potential temperature, 

30/at x 1CT deg K/sec, at 26/1200Z May 1973. 
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900 x 10 gta kgm” /aoc, 26/12Z 850 8w/3t x 10 gta kgta“ l /»cc, 26/222 

Figure 13. The local change of mixing ratio, ciw/dt 
x 10 + gm kgm" /sec, at 26/1200Z May 1973. 
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DISTANCE SECTION A-A' 


The local change of mixing ratio, Qw/gt 
x 10^ gm tcgnT^/sec, for the cross-section 
A-A', 26/1200Z May 1973. 


DISTANCE SECTION A-A’ 


The local change of potential temperature, 
30 /at x lcA deg K/sec, for the cross-section 
A-A 1 , 26/1200Z May 1973. 
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Figure 16. The local change of potential temperature, 

30/3t x 10^ deg K/sec, at the advectlve fore 
cast time 26/152DZ May 1973. 
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Figure 17. The local change of mixing ratio, Sw/iJt 
x 10 4 gm Icgra'Vsec, at the advective 
forecast time 2&/1520Z May 1973. 
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Figure 18. The local change of mixing ratio, $w/dt x 
lCr 1 ' gm kgra" Vsec, for the croee-section 
A-A', at the advective forecast time 26/1520Z 
May 1973. 



Figure 19. The local change of potential temperature, 

B 0/ a t x 10^ deg K/sec, for the cross-section 
A-A', at the advective forecast time 26/1520Z 
May 1973. 
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Figure 20. 


Possible mechanism of squall line develop- 
ment (from Miller, 1972). 



CHAPTER VII 


SUMMARY AND RECOMMENDATIONS 

Several features of the advection analysis with variational 
optimization analysis techniques are especially noteworthy. The first 
is inherent with the variational technique; the governing equations 
and conditions imposed on a model are simultaneously applied. For 
example, the data filtering is accomplished simultaneously in all 
dimensions with the ability to specify the amount of scale smoothing 
for each dimension. 

The method of determination of the vertical component of motion 
is another major feature. The use of the non-homogeneous , anisotropic 
exponential weight function for the objective analysis of the data 
fields produces physically realistic divergence values at all pressure 
levels. These divergence values are then used in a variational tech- 
nique to determine the vertical motion, within the boundary conditions 
imposed, consistent with mass continuity. The vertical motion field, 
crucial to the three Bpatial dimensional analysis of advection of tem- 
perature and mixing ratio, can be assumed to be accurate with a high 
degree of confidence. 

Several promising areas for continued research constitute the 
final major feature of this research. Research should continue to con- 
firm the existence of the vertical circulation pattern identified with 
the 26/1200Z May 1973 case study. This vertical circulation may produce 
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vital clues to the development and maintenance of the severe thunder- 
storm and its parent squall line. 

Another promising area for future research would be the adapta- 
tion of the model to a finer mesh grid over a smaller geographical area 
in order to study the advection fields on the scale of the thunderstorm 
and its parent squall line. The densely- instrumented Beta network, 
operated by the National Severe Storms Laboratory, in central Oklahoma 
is a distinct possibility as the data source. Because the degree of 
scale smoothing is controllable at each analysis step from the initial 
objective analysis to the final four dimensional advective analysis, the 
mesoscale features of interest can be preserved. 

The final recommendation for further research is to investigate 
the short range forecast capabilities of the variational models described, 
In addition to using the advective patterns and forecast fields to iden- 
tify the primary areas of squall line initiation several hours in advance, 
the distinct possibility exists of parameterizing the data to provide 
severity criteria, such as stability indices. These operational possi- 
bilities should be explored in an effort to increase the accuracy of the 
severe weather watches and warnings. 
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APPENDIX A 


COMPARISON OF CHESSMAN AND ANISOTROPIC EXPONENTIAL 
WEIGHTING FUNCTIONS— A CASE STUDY 


Both the Cressman Weight function and the anisotropic weight- 
ing function assign values to grid points from irregularly spaced ob- 
servation points by using the equation 


$ 


ij 


SH k *k 

E % 


where and are the variable values at the observation point and 
grid point respectively and W k is the weight for that observation point. 
Each observation is weighted according to a prespecified relationship 
and the sum of these weighted observations is normalized by the sum of 
the weights. 

Since the Cressraan weighting function is distance dependent 
only, it is isotropic, and because the weighting function is not a func- 
tion of location, it is also homogeneous. The analysis of a wind field 
or any other anisotropic, non-homogeneous field with an isotropic, homo- 
geneous weight function is likely to introduce errors into the objective 
analysis . 

In the early stages of this study, the Cressman weighting func- 
tion was used to objectively analyze the wind fields at each of 19 • 
levels at 50 millibar increments from 1000 to 100 millibars over the 


15 by 21 horizontal grid shown in Figure Al. However, at certain levels, 
the divergence of the horizontal wind was not consistent physically with 
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Figure Al. The horizontal grid established for the study. The grid spacing is 

approximately 125 km. The darkened station circles indicate locations 
of synoptic upper air reporting stations. 
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the synoptic situation. Since the accurate computation of the diver- 
gence is necessary in this investigation, significant error in the 
initial objective analysis could easily produce unknown error terms 
in the remainder of the calculations. 

/ 

As an example, the divergence field was calculated from the 
horizontal wind components assigned by a successive correction objec- 
tive analysis with the Cressman weighting function. The observed data 
field is shown in Figure A2. 

As shown in Figure A2, a closed low height center was in central 
Texas with a sharp ridge extending from New Mexico through Colorado and 
Kansas, Wind maxima were present in the northwest and southeast por- 
tions of the grid. Synoptically , convergence is associated with the 
confluence of airflow and divergence 1 b associated with diffluence. 

The divergence calculated from the Cressman analysis is shown 
in Figure A3. The calculated divergence values in northern Arizona and 
northern New Mexico do not appear to be realistic— apparently the scheme 
treated the wind at INW (Winslow, AZ) and the wind at ABZ (Albuquerque, 
NM) to be strongly convergent. Another apparent error exists in north- 
ern Texas and central Oklahoma where the very high divergence values 
appear to be the result of the wind differences at AMA (Amarillo, TX) , 

FTW (Fort Worth, TX) , and OKC (Oklahoma City, OK). A streamline analysis 
suggests significance error. 

The same field was then analyzed using an anisotropic exponen- 
tial weighting function, with only one pass through the data field. 

The divergence values calculated from these analyses are shown in Figures 
A4a, A4b. and A4c, and appear much more realistic. Notice that the 
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Observed wind field (knots) at 200 mb, 11/1200Z June 1973. The wind field 
indicates wind maxima in the upper left and lower right portions of the grid, 
a closed low height center in the lower center and an inclined ridge from the 
lower left to the upper center. 



63 



ORIGINAL* PAGE IS 
OF POOR QUALUr 


L200Z June 1973. 
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selection of different sets of K and (3 change the magnitudes but not 
the patterns of divergence, 

The divergence values calculated from the Cressman analysis 
have a greater magnitude, even though the patterns were derived from 
wind fields subjected to greater implicit smoothing than in the analy- 
sis based on the new weighting function. Also, a comparison of the 
divergence patterns shows significant differences in sign of the 
divergence as well as the magnitude. The differences are most pro- 
nounced in the vicinity of the smaller wavelength features of primary 
interest to the advection investigations. 

The differences appear to be a function of the isotropy of the 
Cressman weight function; perhaps any isotropic weighting function will 
introduce significant error, which would affect the validity of the 
computations and conclusions of later parts of the study. In any event, 
the necessity of accuracy in the initial objective analysis of the data 
field is emphasized. 



APPENDIX B 


THE USE OF VARIATIONAL OPTIMIZATION ANALYSIS IN DATA FILTERING 

Variational calculus may be used to find the stationary points 
of a functional, an integral of a function which may be composed of 
several dependent variables and their first and second order time and 
space derivatives, in a manner analogous to the use of differential 
calculus to find the stationary points of a function with regard to 
each of its variables. When the functional is in quadratic form, the 
stationary values of the functional will be minimum (Sasaki, 1970a), 

The functional is composed of a series of weighted terms rep- 
resenting conditions to be satisfied, either exactly with strong con- 
straint weighting or approximately with weak constraint weighting. The 
strong constraint weights are termed Lagrange multipliers and are 
treated as introduced dependent variables with respect to the functional. 
The terms representing conditions to be exactly satisfied are premulti- 
plied by a Lagrange multiplier for inclusion in the functional (Sasaki, 
1970a). For terms representing conditions to be approximately satis- 
fied, the square of the term Is multiplied by a variable independent 
of the functional. This term is known as a weak constraint; terms with 
a Lagrange multiplier are known as strong constraints. These weak and 
strong constraints may represent the model dynamics, kinematics, thermo- 
dynamics, and energetics. One weak constraint term common to the 
majority of all functionals is the observational constraint, a weak 
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constraint, which relates the observed or initial field (the "tilde" 
field) to the true field, or that which satisfies all constraints in 
the formalism. By specifying the observational constraint weight, 
the amount of departure from the initial field can be controlled. 

'Ibis ability of the variational optimization analysis method 
to simultaneously impose dynamic, statistical, kinematic, and other 
conditions in both diagnostic and prognostic forms is perhaps the 
major advantage of the variational method (Sasaki, 1971). The most 
important uses of variational optimization analysis have been in data 
filtering and in the initial assignment of data fields consistent with 
the dynamics or other conditions of the model incorporating the data. 

Thus, a geos trophic forecast model can be initialized with wind and 
geopotential height fields adjusted to geostrophic balance. 

Each term or squared expression in the functional is minimized; 
as many terms as can be adequately treated analytically and computationally 
may be added to the functional. However, due to noise and unresolvable 
wavelengths in the data, not all conditions can be simultaneously exactly 
satisfied. The degree to which a condition must be satisfied can be 
specified through the assignment of the weights to the weak constraints. 

A low value of a weight will allow more deviation from the constraint 
than a high value. 

The weak constraint weights are specified independent of the 
functional but may vary as a function of the spatial and temporal coor- 
dinates. The value of the weights may be determined empirically from 
examination of the final or true data fields, from statistical descrip- 
tions of the data fields, from a priori knowledge of the reliability 
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or accuracy or the data, which may be a function of source, time, 
location, or observation method; and from analytic solutions of the 
analysis equations derived from the variational formalism. For the 
simpler finite difference approximations and assumed wave forms of 
the solution data fields, the response function for varying weights 
can be determined. 

In a low pass filter, the desired result is to suppress the 
amplitudes in the higher frequencies in the data field by locally 
altering the values of the data. Tills can be accomplished with varia- 
tional optimization analysis by including first derivative terms as 
weak constraints in the functional (second derivative terms, represent- 
ing curvature, may also be used (Wagner, 1971)). A high frequency in 
the observed data field would appear as locuily large slope values; 
the effect of the weak constraint as the first derivative is to strongly 
reduce locally large slope values and to have little effect on the low 
slope values representing low frequencies in the data, 

The variational functional selected as the filter contains an 
observational constraint and the first derivative constraints in three 
dimensions. The functional is 


J 


fl 


j . ~ 2 

t, ffi(u-u) 



2 + O 1 2 > + “p t ( i> 2 + C> : 2 d * d * « 


where a, o/^ and <y are the weak constraint weights; u and v are the 
analyzed fields and u" and v are the observed or assigned fields. 

The stationary value of the functional is found by equating 


the first variation of the functional to zero. Some basic properties 
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of the variational operator, 6, are: 

1. Variation and differentiation are permutable (SaBaki and 
Lewis, 1970; Sheets, 1973). 

2. Variation and integration are permutable (Sasaki and 
Lewis, 1970; Sheets, 1973). 

3. In finite difference form summed over the domain, the 
commutation property of the variation is (from Sasaki, 

1969; Wagner, 1971): 

4. In the following derivations, each variable $ varies by 
an arbitrary non-zero amount 6$ in the interior, 6$ is 
zero on the boundaries os the boundary condition. 

5. The variational operator acts on all variables in much 
the same way as a partial differential operator. The 
laws of variation of sums, products, ratios, and powers 
are completely analogous to the corresponding laws of 
differentiation (Hildebrand, 1965). 

The derivation of the Euler-Lagrange equation in the finite difference 

form for numerical solution is, where X, Y, and P are the summation domains 

6J - 0 = S 6[a{(u-u) 2 + (v-v) 2 } + ^{(v u) 2 + (nu) 2 + 

XY P r x y 

(V x v ) 2 + (V y v) 2 } + cv p {(v p u) 2 + (v p v) 2 }] (2) 


6J = 0 = ££ 2 2a (u-u)6u + 2 a (v-v) 6v + 0^(2 7 x uV^6u + 


2 V y uv y 6u} + o- f [2 V x vv x 6v + 27 y vv y 6v] 


+ tf p (2 y p u Vp6u + 2 VpV v p 6v} (3) 

Use of the commutation properties of the variational operator produces 

6j = 0 = SEE [aku-u) - a^V^u - a^u - olV 2 u3 26u 
jjy p XX z y PP 

+ fa (v-v) - q - e v 2 v - - a v 2 v) 26v 

L ' f x f y P P 


( 4 ) 
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Since the variation is arbitrary on the interior, the only way 
that 6J can be zero is for the coefficients of 6u and 6v to be zero. 
This condition forms the Euler-Lagrange equations 

«f + ?y U > + °p V p U - au B 'QU (5) 

a f (v^v + yyV) + OpVpV - av - - ofv (6) 


These decoupled equations form the analysis equations. They 
are elliptic partial differential equations which can be solved by 
Richardson or Liebmann relaxation numerical techniques. 

The relaxation technique for specified boundary values is a con- 
vergent, iterative process in the grid interior. The initial fields form 
the first guess field, a residual is computed and a correction is applied 
to the current value, and the iteration continues until there is conver- 
gence to a pre-selected tolerance. In the use of the variational tech- 
nique to filter data fields, the author noted that convergence was 
significantly increased by eliminating sharp gradients on the boundaries 
with a Hanning filter (Blackman and Tukey, 1958). The boundaries then 
remain constant during the relaxation. 

To solve the analysis equations, equations (5) and (6) are ex- 
pressed in center-space finite difference form and then non-dimens ionalized 
and scaled so that the terms are of the same magnitude in the analysis 
equations. The center space finite difference form for a three- 
dimensional grid with the i, j, and k directions correspond to x,y, 


and p respectively. Equation (5) becomes 


d 2 (u i,j-fl,k + U i,j-l,k + U i+l,j,k + U i-l,j,k ' 4u i,j,k 
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where d is the grid distance in the i and J directions and h is the 
grid distance in the k direction. 

Collect terms and rewrite as 

7 < u i,j + i,k + u i, j-i.k + u i+i,j,k + “l-i.i.k* + 

5> 


„2 < u i,j,k+l - “ A 


4a, 

(-a - — 2 


_E v u 
.2 > u i,j,k 


B u 




( 8 ) 


a u 


t.J.k 

thus, Bu. ,+A+C=Ois the desired numerical solution. For the 
n c ^ iterative pass (and notationally dropping the subscripts of u) , 

B u n + A n + C * Residual 11 (9) 


when a correction is applied to u, the residual vanishes to produce 

Bu ftfl +A n +C = 0 (10) 


If equation (10) is subtracted from (9), the result is 


n4-l n Residual 
u = u 


n 


B 

Substitution of B back into the expression produces 

n+1 n . Residual 0 
u = u + 


2 + ^ + 

d 2 h 2 


( 11 ) 


( 12 ) 


The substitution of (9) into (12) produces the correction to apply u 
sequentially at each grid point until the desired convergence is 


n+1 


n , 
u + 


A n + Bu n + C 

a 2 h 2 


achieved. This is 


(13) 
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. X S X \ x*l 

or (5 + — 5 - + — £) u , “ c*u 

d h 1 *J» K 


i,JA + ^2 (u i,j,k + “l.J.k-l* 


+ ,2 (u i,j+l,k + u i,j-l,k + u i+l,J,k + u i-l,J,k> 


(14a) 


« + ■ ? + ^ + + 


and 




+ d 2 (v i,j+l,k + v i,J-l,k * V l+1, j,k + v i-l,j,k ) 


(14b) 


To examine the analytical response function, which is defined 
as the ratio of the analyzed wave amplitude to the initial wave ampli- 
tude as a function of wavelength, assume that u and u have the analytic 
form 


~ 7 i(kx + ky + np) 

u *» A e 

(15) 

u e A e i(kx + ky + np) 

(16) 


where k and n are wave numbers in the horizontal and vertical directions, 
and substitute Into the analysis equations (5) and ( 6 ). For the 
analysis equation for u only, after the substitution of (15) and (16) into 
(5), carrying out the indicated differentiation (in continuous form), 
and canceling the exponential term which remains unchanged, the result 
is 

- 2k 2 Q' f A - n 2 a p A - c? (A- A) (17) 

Expressing the analytic response function, R, as A/A, the response 

function is 22 

2k a* n ql 1 

R - (1 + “ “ + -7T 2 ) (18) 

a a 

To determine the filter response function of the finite differ- 
ence analog, we need initially consider only the response in the x 
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direction. The assumed forms of the analyzed and initialized u fields 
are, in finite difference form 

u - A e 1<krAlt> (19a) 

u - Ae 1(ltriix) (19b) 

where k “ 2tt/L, the wave number, and x ■ rAx, and Ax « grid spacing. 
The analysis equation (7) in the x direction is 

~2 ^ u i+l + u i-l " 2u i^ " ( 20 > 

Ax 

Substitution of equations (19a) and (19b) into (20) produces 

— r [exp(ik(nfl) Ax) + exp(ik(r-l)Ax)-2 exp(ikrAx)] 

Ax 

*» O'(A-'K) exp(ikrAx) (21) 

Cancellation of the exponential term yields 

(e ik6x + e ik & - 2) -£(A-X) (22) 

Ax 


Use of the identities 

2 cos(x) = e*"* + e ix 

2 cos (2x) - 2 «= -4 sin 2 (x) 


produces 


~ 2 (-4 sin 2 (kAx/2) ) 
Ax 


a 

A 


(A-A) 


which simplifies to the response function 

A „ , .2 


R 


A j 2 . tt . . -1 

_ 3 (1 + — ~ 2 S i n ( - Ax)) 

A aAx 


(23) 


(24) 
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Similarly, the response function in the p plane is 

R “ + “T sin2 < ? ^p ^' 1 (25) 

P »AP p 

The low pass filter characteristics of the first derivative 
filter may easily be verified from equation (24). As the wavelength 
L, expressed as a multiple of Ax, increases from twice the grid spacing 
(which corresponds to the Nyquist frequency) to large values, the 
filter response function approaches a value of lj i.e., no amplitude 
change for that wavelength, with the majority of filtering tech- 
niques, aliasing will occur for frequencies greater than the Nyquist 
frequency, the highest resolvable frequency. Since the Nyquist fre- 
quency is defined in terms of the grid spacing, the filtering charac- 
teristics are a function of the grid spacing. is advantageous 
therefore to select a grid interval smaller than the data measurement 
point spacing. This results in a greater discriminatory filtering 
capability . 

To this point, the equations and derivations may be regarded 
as either dimensional or as scaled and dimensionless. However, the 
computations and the selection of the weights are simplified when the 
equations are in scaled, non-dimensional form. 

For the analysis equation (5) 

.2 2 . . 2 . ~ «■“> rs* 

a f <V + Vy u -> + Op Wp u ) * a U = ■ Of u 

let each variable be the product of a "star" term and a "prime" term 
where the primed term is dimensionless and of the magnitude of one. 
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Thus, 


and, 
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= 
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The substitution of the variables 
the collection of all dimensional 


a 2 1 2 1 

(v K -• ♦ v* «■) 


and operators into equation (5) ar.d 

terms produces 
* 


o' 


+ > 


-~*.2' 
a h 


Vp 


u' = y' (u'-u') (26) 


The selection of the scale values is largely independent on the 
scales of interest in the study and the grid spacing established. For 
this study, the scale values are 

u* = 10 m/sec 
d = 125 hr. 
h e 50 mb 


To reduce the left hand side of equation (26) to the order of 



1 

d 


2 


h 


2 


magnitude of one, select 
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Equation (26) then becomes identical to the form used in 
derivation of the filter response. The only difference is that vari- 
ables and parameters have scaled, non-dimensional values. For example, 

2 2 

in equation (24) , the dimensional value of Ax is near (125 km) 5 the 

2 2 

non-dimensional value of Ax is (1) . 

The response function of equations (24) and (25) may easily be 
calculated for any wave length in multiples of the grid spacing. The 
response functions for both first derivative filters and second deriva- 
tive filters are published by Wagner (1971). The response function for 
the first derivative filter term comprises a family of curves. Since 
the measurement point spacing (not the grid spacing) in the horizontal 
and vertical defines the minimum resolvable wavelengths in the data, 
the filtering response should differ for horizontal and vertical wave- 
lengths. This is controllable by adjusting the ratio of the non- 
dimensional weights, as shown in Wagner's study. The response functions 
shown in Figure Bl illustrate the family of response functions for this 
inves tiga tion. 

It should 1: e noted that this filter technique does not produce 
a phase shift or frequency change. This feature is found in other 
finite difference filters (Haltiner, 1971) . 


RESPONSE 
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Figure Bl. Response function of the first derivative filter 

as a function of wavelength, expressed as multiples 
of the grid interval, for varying ratios of filter 
weight, o', to observational weight, a?. 


APPENDIX C 


SEVERE THUNDERSTORM DEVELOPMENT THROUGH THE WARM, DRY INTRUSION 

The initial conditions at 26/1200Z May 1973 and the advective 
forecasts from the variational model indicated that at the time and 
location of squall line development, the low-level moist air was over- 
lain by a deep layer of warm, very dry air. The presence of the 
thunderstorms thus indicates an indirect circulation; cooler moist 
air is rising convectively through the warmer, dry layer above. At 
cloud base, the thunderstorm updraft air would be negatively buoyant 
with respect to the mesoscale environmental conditions. The existence 
of this indirect circulation is supported by an examination of the 
ATS III imagery for 26/1200Z May 1973 in conjunction with the varia- 
tional model results. Six selected frames with brief descriptions 
are included as Figures Cl to C6. 

Figure Cl depicts the pre- thunderstorm environment. The warm, 
dry air of the dry intrusion extends eastward over the stratus in 
eastern Oklahoma (apparent also in Figure 3, the Tinker AFB sounding) 
from the surface position to the west of the stratus boundary. The 
initial thunderstorm development occurred on the western boundary of 
the moist air. This is a preferred location for squall line thunder- 
storm development (Tegtmeier, 1974), which agrees well with Schaefer's 
model of the relation of the surface dryline to the inversion marking 
the lower boundary of the dry intrusion (Schaefer, 1974), 
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As the thunderstorms develop, the warm, dry air continues to 
advect toward the east, as shown in Figures C3 to C5. This eastward 
advection is indicated also in the forecasts from the variational 
advective model as depicted in Figures 15 and 19, for example. 

The mesoscale environment of the severe thunderstorms shown 
in Figures C5 and C6 appears to be a low level moist layer (below 850 mb), 
capped by an air inversion marking the boundary between the moist air 
and the much warmer dry air intrusion, which extends at least through 
the 700 mb level. The mesoscale environment thus indicates an indirect 
convective circulation below the 700 mb level; cool moist air rises 
convectively through warmer, dry air. 

As the thunderstorms develop to severe levels, they signifi- 
cantly modify their mesoscale environment in the vicinity of the storm 
and in the area downwind from the storm, since the storms typically 
translate at a rate less than that of the environmental winds. The 
scale interactions which exist at the levels of the dry air intrusion 
immediately upwind of the storm may be a significant factor in the 
development of the storm to severe levels. The variational advective 
forecast results Indicate that the intrusion is very dry, with virtually 
an unlimited capacity for evaporation, and warm with a near neutral 
lapse rate. Thus, evaporation of cloud material into the dry air would 
produce a strong subsistence of the cooler air, especially since even 
with a dry adiabatic descent, the equilibrium level may be several 
hundred millibars below the original level of cooling. As surface 
heating warms the cooler layer beneath the inversion, this subsiding 
air may reach the surface. 

A major effect of the warm, dry intrusion in he severe storm 
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environment is to facilitate a detrainmcnt of the storm air into the 
mesoscale environmental air, as opposed to the entrainment of environ- 
mental air into the storm. The zone where rapid evaporative cooling 
occurs will protect the storm updrafts from the dry environmental air, 
and perhaps is necessary for development to severe levelB. 

Recent published research on severe storms support detrainmcnt 
into the dry intrusion (Davies -Jones , 1974; Marwitz, 1972a and 1972b; 
Foote and Fankhauser, 1973) . Davies-Jones observed that soundings in 
severe storm updrafts indicated that the cores of strong updrafts were 
undiluted by the environmental air. He also observed that the sounding 
balloons rose almost vertically in the storm in spite of pronounced 
environmental shear. A negative buoyancy at cloud base was observed 
from aircraft measurements by Marwitz. Based on time-section analyses 
of two storms, Marwitz inferred that the descending air to the south- 
west (the windward side) of the thunderstorms had subsided from above 
the cloud base. He suggested that this mid- tropospheric air is directed 
downward as a result of evaporative cooling and deflection; this sub- 
siding air reaching the surface layers could intensify the gust front, 

Foote and Fankhauser (1973) state that negatively buoyant in- 
flow air is a fairly general observation for large, persistent storms 
and conclude that dynamic forcing drives the low-level inflow into the 
storm. Correlations of the surface conditions, as indicative of the 
low-level inflow, with squall line thunderstorm development have been, 
good. Sasaki (1973) concluded that intense surface moisture conver- 
gence coincides to squall line development and local severe weather 
areas and that a local downward momentum transport is very important 
to squall line formation. 
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The role of the dry air intrusion to severe Btorm development 
is as yet speculative. The scale interactions are difficult to param- 
eterize, describe, and model and the measurements made of severe thunder- 
storms may sample either the mesoscale unaltered environmental conditions, 
the thunderstorm scale conditions, or the interacting zones. Some con- 
clusions may be hypothesized, however. The warm, dry intrusion facili- 
tates detrainment with a resulting downward transport of modified air 
and deflection of environmental air. This mechanism provides for a 
protection of the storm interior from the dry air and the environmental 
shear through the dry layer. The subsiding air provides for a downward 
momentum transport and will serve to intensify the gust front. Another 
significant factor not directly commented upon is that the inversion 
layer marking the lower boundary of the dry intrusion serves to cap the 
potential instability of the moist air below. The relation of the dry 
intrusion to the initial development of thunderstorms was not indicated 
in the variational advective model results. 




Figure Cl. ATS III, 26/1819Z May 1973. The location of th» dry 
air intrusion is approximated by the clear area, in- 
dicated by A. Ihe cloudiness IJ is associated with 
the southeastward moving cold front from the low 
pressure center in central Kansas. H»e cloudiness 
C is stratiform beneath the dry air intrusion. There 
is no thunderstorm activity in Oklahoma at this time. 



Figure C2 . ATS III, 26/1858Z May 1973. The location of the 

first reported thunderstorms in Oklahoma is indicated 
by D. 
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Figure C3. ATS III, 26/2006/. Mav 1973. The thunderstorm 

activity in Oklahoma has developed southward and 
is aligned as the well-defined squall line E. 



Figure C4. 
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ATS III, 26/2031Z May 1973. The rapid development 
of the squall line E in the 25 minutes between 
Figure C3 and Figure C4 is very apparent. The first 
of the numerous tornadoes with this squall line was 
reported at 2052Z with the thunderstorm in Kansas 
Just north of the Oklahoma border. Tornadoes were 
also reported at 2051/ with the thunderstorm complex 
C in southern Kansas. 
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Figure C5. ATS III, 26/2110Z May 1973. The squall line E is 

now producing funnel clouds and tornadoes alcng its 
entire length. Moderate thunderstorms are occurring 
along the cold front in northwestern Oklahoma and 
the thunderstorm complex in central Kansas. 



Figure C6. ATS 111, 26/2201Z May 1973. The fatality-producing 
tornudo occurred near 22.00Z at F in the southern end 
of the squall line. Large tornadoes were also reported 
near 2200Z at G. 
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